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ABSTRACT 

We perform fully-kinetic particle-in-cell simulations of an hot plasma that expands radially in a 
cylindrical geometry. The aim of the paper is to study the consequent development of the electron 
temperature anisotropy in an expanding plasma flow as found in a collisionless stellar wind. Kinetic 
plasma theory and simulations have shown that the electron temperature anisotropy is controlled by 
fluctuations driven by electromagnetic kinetic instabilities. In this study the temperature anisotropy is 
driven self-consistently by the expansion. While the expansion favors an increase of parallel anisotropy 
(Tii > T±), the onset of the firehose instability will tend to decrease it. We show the results for a 
supersonic, subsonic, and static expansion flows, and suggest possible applications of the results for 
the solar wind and other stellar winds. 
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1. INTRODUCTION 

The steady-state solution of the most simple dynam- 
ical equations that describe the evolution of a stellar 
wind dates back to the hydrodynami c model for th e so- 
lar wind proposed by Parker in 1958 (Parker 1958}). In- 
creasingly sophisticated models have been developed to 
take into account the different observed characteristics 
of the solar wind, and modeling has been based on ana- 
lytic, semi-analytic and simulational methods \Hollwea\ 
2008). The magnetohydrodynamics (MHD) approach 
is widely used for global modeling, but since it treats 
the plasma as a fluid it does not include any effects due 
to non-Maxwellian particle velocity distributions. More 
recently, kinetic simulations have tried to take in ac- 
count and explain features related to the observed non- 
thermal distribution of parti cl es in velocity spac e (see e .g. 
\Landi and Pantellini (|2003ft : \Zouaanelis et al\ (|2005l ) ). 
The aim of such exospheric models is to provide a real- 
istic description of the wind dynamics that includes the 
transition from a collision-dominated to a collisionless re- 
gion. In doing so, however these models do not include 
the effect of plasma instabilities and therefore cannot be 
regarded as completely self-consistent. 

For the solar wind, the importance of small scale fluc- 
tuations, associated with kinetic plasma instabilities gen- 
erated by non-Maxwellian particle distributions, is now 
widely recognized. This has come about through a con- 
vergence of observations, theory and simulations. It is 
argued that many macroscopic quantities that charac- 
terize the solar wind, such as the particle temperature 
anisotropy or the electron heat flux, are always observed 
with values that are bounded by the possible onset of 
kinetic instabilities. The general argument is the follow- 
ing: the expansion of the flow leads to a distortion of the 
distribution function, which represents, for example, an 
increase in the temperature anisotropy. The distortion of 
the distribution function can be enough to trigger linear 
instability, i.e., there is some free-energy available which 
can create plasma fluctuations. Based on simulations of 
the initial value problem for such instabilities, it is evi- 
dent that the fluctuations act to reduce the distortion of 
the distribution function. In other words, the primary 
effect of the fluctuations produced in the instability is to 



restore a stable (or marginally stable) distribution func- 
tion. In the case of a temperature anisotropy driven by 
expansion, we might expect the onset of instabilities as- 
sociated with temperature anisotropy to act as an upper 
limit for the anisotropy. This argument can be broadly 
applied to many situations, such as formation of heat 
flux, or indeed compression flows. 

In the solar wind context it is expected that the ex- 
pansion of the wind away from the Sun produces a par- 
allel temperature anisotropy T» > T± (with reference 
to the magnetic field direction) due to the conservation 
of adiabatic invariants. If one adopts a flui d viewpoint 
(for i nstance the 'double-adiabatic' theory \Chew et ali 
119561 ) that is often used in this context), the anisotropy 
is not bounded and it should increase indefinitely for in- 
creasing distance from the Sun. However it has been 
observed that the tem perature anisotropy is limited 
in value both for ions (I Gary et aZJ 1200 J: \Kasver et all 



120021 : \Hellinaer efaH 120061: \Matteim let all 12007ft and 
electrons ( Gary et al\ l2005t \Stverak et ali [2008ft . and 
it has been suggested that this is due to the onset 
of kinetic instabilities, such as the firehose for Tji > 
T± and the mirror or ion-cyclotron instabilities for 
Tii < Tj_. Recently, several computer simulations have 
been able to show that the kinetic instabilities are ef- 
fectively able to constrain the growth of temperature 
anisotropy. Those inc l ude hybrid simulatio ns for protons 
[Hetti nger et a l. 2003; Matteini et at 200 6]) and particle- 
in-cell (PIC) simul ations fo r electrons ( Gary and Wan 
1996t [Gary et all [2000|: \Gary and Nishimyral 1200 
Camporeale and Buraesd 2008ft . 



Most simulations study the initial value problem for 
an instability, starting with some initially unstable 
anisotropy value and following the evolution of the sys- 
tem to marginal stability. These simulations are usually 
interpreted in the framework of the linear theory of the 
Vlasov-Maxwell equations. The linear theory is devel- 
oped for a non-expanding plasma embedded in a uniform 
magnetic field, in a periodic Cartesian geometry. It is 
this theory which is used to determine the marginal sta- 
bility boundary used to confirm the role of linear insta- 
bilities as constraining the observed particle anisotropy. 
In other words, the paucity of observed periods of large 
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anisotropy is interpreted as a consequence of fluctuations 
driven by anisotropy instabilities which stop the plasma 
reaching a high anisotropy state, or rapidly relaxing it 
if it does. The bounds of the anisotropy are calculated 
using a non-expanding Cartesian geometry. 

All the simulations performed so far do not include 
the radial expansion of the wind in a completely consis- 
tent manner. Because of the computational difficulties of 
simulating an expanding plasma in a fixed frame of ref- 
erence, simulations using an 'expanding box model' have 
been developed. In this type of simulation the idea is to 
use a Cartesian computational domain, but to 'stretch 
out' the domain as the plasma expands, thus distorting 
the simulation box in at least one dimension. As a con- 
sequence, the equations of motion have to be modified, 
taking into account the inertial forces due to the expan- 
sion of the box, and the coordinates must be continuously 
rescale d. For plasma s i mulat ions this method was initi- 
ated in I Gravvin et all (|1993|) for MHD, and it has been 
successively applied to hyb rid simulations (with fluid 
electrons and particle ions) m \Hellinaer et cd\ (|2003l ). Al- 
though the simulation results have been relatively suc- 
cessful, and shown to be con sistent with satellite obser- 
vations \Matteini et aLll2007t ). the expanding box model 
is still unable to treat the expansion consistently. It can 
be argued that the rescaling of the box effectively pro- 
duces an unphysical mode coupling due to the allowed 
wave vectors changing continuously over time. Energy 
that resides in certain modes at one time is artificially 
channelled towards other modes, as time evolves. The 
counter argument is that this mode evolution is actually 
expected as precisely a result of the expansion. 
In this paper we present the first fully-kinetic PIC sim- 
ulations, with realistic proton-electron mass ratio, of a 
plasma that expands radially in a decreasing magnetic 
field, in a fixed frame of reference. The expansion is ra- 
dial in a cylindrical geometry. The scope of this work is 
to study and quantify in a more consistent way the com- 
petition between the growth of parallel anisotropy due to 
the expansion, and the possible onset of kinetic instabili- 
ties. The simulations presented in this paper are 2D-3V, 
i.e. two dimensional in space and three dimensional in 
velocity space. The magnetic field is also 2D, with axial 
component zero. Clearly, this is just a first step towards 
more challenging, and realistic 3D simulations. But, we 
will show that even with the use of cylindrical geometry 
it is already possible to capture the increase of anisotropy 
and the subsequent development of kinetic instabilities. 
With the scales currently feasible, our results are relevant 
to the evolution of the electron particle distribution. We 
will show and compare three different cases: static, sub- 
sonic, and supersonic flow. 

The paper is organised as follows. The feasibility of run- 
ning the simulations has been greatly enhanced by using 
an implicit scheme, and we discuss the main features of 
the algorithm in Section 2. Section 3 is devoted to de- 
scribe the details of the box geometry and size, and the 
characteristic plasma timescale and lengths. We show 
and discuss the results of different runs in Section 4, and 
draw our conclusions in Section 5. 

2. THE IMPLICIT CODE 

The code used in this work is a fully-kinetic, implicit, 
parallel, electromagnetic PIC code called PARSER. Its 



features are described in detail in Markid is et al\ (2009) , 
but for completeness we give here a brief desc ription of 
the al gor ithm. We refer the reader to , e.g., \Pritchem 
(2003) or \Hockney and E astwood} (|1981| ) for a more gen- 
eral tutorial on PIC techniques. 

The algorithm is implicit in time both for the par- 
ticle mover and for the field solver, and it adopts 
the so-called 'implicit m o ment method', introduced 
by \Brackbill and Forslund\ (119821) . an d succ essively re- 
elabor ated in lVw and BrackbilH (|1992j ) , and \Ricci et~al\ 
(2002). The following equations (in CGS units) for the 
conservation of density and momentum are used to ex- 
trapolate the values of the charge density p, and the cur- 
rent density J at a new timestep: 
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where superscripts indicate timestep, and all other sym- 
bols are standard. The closure of equations (jTJ^J) is pro- 
vided by approximating the divergence of the pressure 
tensor at time i + 1 with the value at time i, that is 
V • P i+1 ~ V • P*. From Eqs. (HE]) one can formu- 
late p z+1 and J I+1 as functions of the electric field E l+1 . 
By using those relations in Maxwell equations, and after 
some algebra one ends up with a linear equation for E l+1 
as a function of only old quantities: 

(c6At) 2 [-V 2 E l+1 - VV • (m' ; E 4+1 )] + (I + //)E 4+1 = E ,: 



where the following terms are defined: 
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and the subscript s indicates the species. 
Equation ([3]) is solved by a matrix-free General- 
ized Minimal Residual (GMRcs) iterative linear solver 
$Saad and Schultz 1986). Once the electric field is known 
at time the magnetic field B is advanced using Fara- 
day's law: 

B i+1 = B' — cAtV X E* + *. (8) 

The particle mover pushes the particles to a new position 
and velocity, according to the equations: 
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Equations (|9I10[) are also solved iteratively with a 
Predictor-Corrector technique. Finally, to ensure that 
the continuity equation is satisfied, the electric field must 
be corrected with: 



E. 



E, 



Id 



, V 2 



V ■ E 



old 



Airp. (11) 



3 



The main advantage of the implicit method over an ex- 
plicit one is that it enables a choice of timestep At and a 
grid size Ax that do not satisfy the Courant stability con- 
dition cAt/Ax < 1. This is instead replaced by the less 
stringent accuracy condition v e At/Ax < 1, where v e is 
the electron thermal velocity. The timestep is still small 
enough to resolve the electron gyromotion. Of course 
this benefit is paid for in terms of the computational com- 
plexity of the algorithm, however for many situations the 
possibility of choosing a fairly large timestep results in a 
positive payoff for the total simulation runtime. 

3. SIMULATION SETUP 

We simulate an ion-electron plasma, with physical 
mass-ratio, i.e., rrii/m e — 1836. The plasma expands 
radially on a 2D disc, and the magnetic field is forced to 
have zero axial component so it is a 2D vector field. The 
geometry of the box is shown in Figure [U The grid is 
Cartesian (x,y) and covers the trapezoid ABCD. The 
oblique sides AB and DC form a 90° angle, and we apply 
periodic boundary conditions on these sides. Therefore, 
a complete plane geometry is recovered by applying three 
successive 90 degrees rotations of the box. This four-fold 
symmetry is used to reduce the computational effort, and 
does not affect the short wave length fluctuations that 
develop. 

To ensure a correct periodicity along the azimuthal di- 
rection, the particles that escape the boundary DC are 
re-injected from the boundary AB, and vice versa. In 
doing so, their trajectory and velocities must be appro- 
priately rotated by 90° . We show an example, in Figure 
[TJ where a particle moving from point 1 to point 2 is re- 
injected at point 4, as if it was coming from point 3. The 
vector velocity (v x , v y ) on the plane must also be rotated. 
The same argument applies to the electric and magnetic 
fields on the boundary, where the x-component on the 
AB side is imposed to match the y-component on the 
DC side, while the y-component on AB must be equal 
to the x-component on AB, with a change of sign (i.e. 
E y on AB is equal to — E x on DC). The out-of-plane z 
component is treated as periodic in the standard way. 

We define three different regions in the box: an inner 
boundary region AEGD, an active region EFHG, and 
an outer boundary region FBCH. Particles are initially 
loaded in all the three regions, with a density that varies 
as 1 jr, with r the distance from the origin (0, 0) (which is 
out of the box) . Except for the static case, they are ini- 
tialised with an isotropic Maxwellian distribution, and 
with a radial mean velocity V m . The initial magnetic 
field on the plane is also radial, pointing outwards, and 
it decreases in magnitude as 1/r. The initial electric field 
is null. While particles are allowed to move anywhere in 
the box ABCD, the field solver advances the fields only 
in the active region. This effectively produces boundary 
conditions on the arcs EG and FH, where any pertur- 
bation of the initial fields is forced to smooth out. 

The fact that we have a reservoir of particles in the 
inner and outer boundary regions that move consistently 
with a static electromagnetic field avoids spurious bound- 
ary effects on the sides EG and FH. We are therefore 
mainly interested in what happens inside the active re- 
gion, which has consistent boundary conditions for both 
particles and fields on all its sides. The true boundaries 
of the computational box are of course the sides AD and 



BC, which are sufficiently far from the active region. 
Here, as we said, the electromagnetic field is static; par- 
ticles are allowed to escape, and they are re-injected at 
every timestep on both sides. The re-injection routine is 
computationally expensive, but it mimics the existence 
of a population of Maxwellian particles outside the box, 
with drift velocity equal to V m parallel to the magnetic 
field, and with density that again goes as 1/r. 

The plasma beta for the species s is defined as usual 
as P s — ^jr?S where n is the density, T the tempera- 
ture, and B the magnetic field. It increases linearly with 
the distance from the origin, and the initial values of 
n, T, B are chosen so that /3 — 7.7 on the arc EG, and 
/3 = 15.4 on FH, for both electrons and ions (that is 
T e = Ti). These relatively high values of beta allow the 
expanding plasma to reach quickly a parallel anisotropy 
large enough to trigger the electron firehose instability. 
According to the linear kinetic theory, the only two pa- 
rameters that control the growth rate of an anisotropy 
instability are the parallel bet a j3\\ and the temperature 
ratio . For the electrons, \ Camvoreale and BuraesA 

( 2008} have found that the relationship 



is valid at the threshold of the instability. 

Velocities are normalized to the light speed c, and the 
initial thermal velocities for electrons and ions are respec- 
tively v e — 8-10~ 3 , and Vi — 1.9 -10~ 4 . The subsonic and 
supersonic case that we discuss in the next section are re- 
spectively for V m /v e = 0.625, and V m /v e = 2.5, while for 
the static case V m = 0, but the initial electron anisotropy 
is T±/Tu = 0.7. The box has a maximum of 540 cells in 
the x direction and 1170 in the y direction. The active 
region has a radial extension of Lu>i/c = 1.25, where Wj is 
the ion plasma frequency. The total box has a maximum 
radial extension of 4.22c/uv The box is therefore suf- 
ficiently large to capture waves with wavenumber kc/cOi 
greater than about 5. The timestep is At — 0.05W" 1 , 
and the cell size is Ax — Ay = O.OObbcoj^ 1 , making the 
Courant parameter cAt/Ax equal to 9; the advantage of 
the implicit scheme is evident. 

4. RESULTS 

In the CGL 'double adiabatic' description of a plasma, 
the quantities -g* and % are constant. It follows that 
the temperature anisotropy varies as Tu/T± ~ n 2 /B 3 . In 
our configuration, where density and magnetic field de- 
crease as 1/r, Tj|/Tj_ will grow linearly with the distance. 
Since we start with an isotropic particle distribution, the 
anisotropy is expected to increase until the local plasma 
parameters are such that Eq. (fl2|) holds. At this point 
the electron firehose instability is triggered. 

The linear dispersion relation o f the electron firehos e 
instability yields two solutions \Li and Hab bal 2000). 
One branch is a propagating slowly growing mode, with 
angle of propagation ranging from 0° to about 70°, and 
the other is a non-propagating fast growing mode, with 
wavevector forming an angle between about 30° and 90° 
with the magnetic field. The latter is generally domi- 
nant, but it has been shown that depending on the an- 
gle of propagation and the level of anisotropy, there are 
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wavevectors fo r which the growth rate of the two modes 
is comparable \Camvoreale and BuraesM 2008) . 

4.1. Supersonic case: V m /v e = 2.5 

In order to evaluate the role of instabilities in reduc- 
ing the electron temperature anisotropy in the expanding 
plasma, we compare in this section two simulations. One 
has the setup described in the previous section, i.e the 
box is divided in three regions, and the electromagnetic 
(EM) field is solved only in the active region. The second 
simulation is identical in all the parameters, except for 
the EM field that is kept static through all the box. In 
this way, the differences between the two runs are clearly 
due to feedback effects caused by self-consistently gener- 
ated electromagnetic fluctuations. We will call the two 
simulations 'self-consistent' and 'test-particle' runs. The 
focus of our interest will be the development of electron 
parallel anisotropy. Firstly, however, we show in Figure 
[H the development at three successive times of the am- 
plitude of the magnetic field fluctuations <5B/B (where 
d~B = |B — Bo |, and Bo is the initial magnetic field), for 
the self-consistent run. These images have been obtained 
by successive rotations of the box. At time Tuii = 18 (top 
panel), there is not yet any large scale structure evident, 
and at that time there are many modes at different ori- 
entation but comparable amplitude present, creating an 
almost random patchwork of magnetic field fluctuations. 
As time evolves (middle and bottom panels), a struc- 
ture of fluctuations aligned with the background mag- 
netic field emerges. This is consistent with the develop- 
ment of quasi-perpendicular waves, that at those times 
have superseded the more parallel modes. By performing 
a Fourier transform in time, we have indeed confirmed 
that the waves observed in Figure [2] are non-propagating 
(or alternatively moving very slowly, compared to the 
total time of the simulation) in the azimuthal direction. 
Movies of the evolution show that the structures do not 
propagate azimuthally, but do have features indicating 
that there is underlying convection outwards in the ra- 
dial direction. 

This is a first evidence that a plasma that in- 
creases its electron temperature anisotropy in a self- 
consistent expansion actually triggers a non-propagating 
firehose instability. This is not a trivial result, be- 
cause although predicted by the linear kinetic the- 
ory in Cartesian geometry, this has never been con- 
firmed before by computer simulations, for an ex- 
panding plasma. As we discussed in the Introduc- 
tion, what was already known is that an instability 
would develop if the plasma was starting with a suf- 
ficien t ly anisotropic temperature iGarv and Nishimura] 
120031 : iCamporeale and Buraes~s\ 120081 ). while here we 
started with an isotropic plasma and let the anisotropy 
grow self-consistently, due to the expansion. The phys- 
ical behaviour here is much more complex than for a 
non-flowing plasma in a constant magnetic field. The 
magnetic fluctuations created by the instability are now 
convected outwards with the flow. This effect helps to 
lower further the temperature anisotropy at larger dis- 
tances, by increasing the particle scattering. 

The top panel of Figure |3] shows the development of 
the anisotropy at three different radial distances, where 
the temperatures are averaged over the azimuthal direc- 
tion. The three solid lines are for r = 1.50 (red), r = 1.86 



(blue), and r = 2.22 (black), for the self-consistent run 
(r is normalized to the ion inertial length coj^ 1 ). The 
corresponding dashed lines show the development of the 
anisotropy at the same radial distances, for the test- 
particle simulation. We have also included, for r = 1.50 
only, the value of anisotropy predicted by the CGL the- 
ory (dot-dashed line). The CGL prediction is in good 
agreement for values of T«/T± less than about 2, but 
then it greatly overestimates the anisotropy (also for the 
test-particle run). The effect of the firehose instability 
in reducing the anisotropy at different distances is clear, 
and it is not surprising that the anisotropy is higher when 
the particles are not scattered by electromagnetic fluctu- 
ations. We show in the bottom panel of Figure [3] the 
ratio between the values of solid lines and dashed lines. 
This represents the reduction of anisotropy due to the 
presence of fluctuations. One can see that this reduction 
stands between 20% and 40% depending on location and 
time. 

We want now to quantify the effect of the instability 
by looking separately at the parallel and perpendicular 
temperatures. These are shown respectively in the top 
and bottom panel of Figure |4l where the temperatures 
have been normalized to their value at time T = 0. The 
colours are the same as for Figure |3l and solid and dashed 
lines are again for self-consistent and test-particle runs, 
respectively. Here we have also pointed, with dotted ver- 
tical lines, the approximate time at which, for different 
locations, the local anisotropy and plasma beta are such 
that Eq. (fT2~|) predicts a marginal stability state. In other 
words, after the time denoted for each curve by the dot- 
ted line, the plasma is unstable at that location. It has 
to be reminded however, that Eq. (|12[) has been derived 
from the linear theory of plasmas in an homogeneous 
field, and in a periodic Cartesian geometry. It appears 
indeed that the curves are not sensitive to the particu- 
lar times designed by the dotted lines. We will give two 
possible explanations later. First we comment on the 
general features that emerge from Figure HI Recall that 
the expansion (without the effect of instabilities) should 
result in a monotonic decrease of Tj_, and a constant 
Tii. This is indeed what happens for the dashed lines. 
At a certain time, the value of Xj_ approach an asymp- 
tote, and this happens for successive times at larger dis- 
tance. This is because the perpendicular temperature 
is inversely proportional to the distance travelled by the 
particles. This distance will increase at the beginning of 
the simulation, until it will reach the maximum possible 
value at each location, equal to the distance from the in- 
ner boundary, and than will stay constant. On the other 
hand, the parallel temperature stays roughly constant, 
in the test-particle run, while the effect of the instability 
is supposed to decrease it. Two comments are in or- 
der. First, the stationary value reached by T± is higher 
in the self-consistent run, meaning that electromagnetic 
fluctuations act to reduce the decrease of perpendicular 
temperature. However, no particular change is appar- 
ent in either Tu or T± when the plasma becomes un- 
stable. Second, the parallel temperature suffers a mild 
decrease in the initial stage and then increases at suc- 
cessive times, at larger distance. This is probably due 
to the concurrent saturation of firehose modes that once 
damping tend to enhance the parallel temperature, and 
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the creation of non-thermal features in the particle dis- 
tribution function, as we will show later. Since those 
modes are convected outwards, the increase in parallel 
fluctuation temperature moves outwards too. 

What emerges therefore is that the electromagnetic 
fluctuations are responsible for slowing down the rate at 
which the parallel anisotropy grows. However, we have 
been unable to unequivocally identify the firehose insta- 
bility as entirely responsible for this behavior. As we 
anticipated there are, in our view, two possible explana- 
tions. One is that the linear theory result summarised 
in Eq. (|12|) is not completely applicable when the ge- 
ometry is radial, and the magnetic field and density are 
not uniform. This is a point that deserves a deeper in- 
vestigation, but it could be that the results developed 
over the years for homogeneous plasma are not straight- 
forwardly applicable to more realistic situations (i.e. for 
non-constant magnetic field and density, and not periodic 
structures in the radial direction). A second possibility 
is that small noisy fluctuations, unrelated to the fire- 
hose instability, could scatter the particles and decrease 
the parallel anisotropy even before reaching a linearly 
unstable condition. Moreover, the increase in parallel 
temperature at later times, seems to be an artificial ef- 
fect due to the development of non-Maxwellian features 
in the particle distribution function, such as high energy 
tails. We show in Figure [5] the contour plot of the elec- 
tron distribution function, averaged in the whole active 
region, at times Tui = 18 (top panel), Tuii = 60 (mid- 
dle), Tu>i — 100 (bottom). One can notice that as the 
perpendicular temperature is reduced, the distribution 
becomes asymmetric in the parallel direction. This is 
reflected in the increased parallel temperature. The dis- 
tribution function for states close to the equilibrium is 
therefore non-bi-Maxwellian, and this is consistent with 
many satellite observations. 

4.2. Static case: V m = 

In this section we show the results of one simula- 
tion, where the particles have no initial mean veloc- 
ity [V m = 0), but the initial anisotropy T±/T« = 0.7. 
In this way the firehose instability is triggered in the 
whole active region from the beginning of the simu- 
lation. Similar simulations, but for an homogeneous 
plasma in a double periodic bo x have been performed in 
\Camporeale and Buraesd (|2008h . The purpose of this run 
is to check that although the geometry is 2D in space, and 
therefore some simplifying assumptions have been made 
on the initial profile of magnetic field and density, and 
on somehow artificial boundary conditions, the simula- 
tions are still able to capture the firehose instability and 
the consequent decrease of anisotropy. Moreover, with 
this run we will be able once again to clearly identify the 
decisive role of the expansion. 

We show in Figure |2] the development of <5B/B . Here 
again the dominance of quasi-perpendicular modes at 
later times is evident. The growth rate of the instabil- 
ity is now higher at larger distances, where the plasma 
beta is higher. Since we start from an unstable plasma, 
and the counter-effect of the expansion is now absent, 
the fluctuations reach amplitudes slightly higher than 
for the supersonic case (Figure (2). If we look at the 
temperature anisotropy (Figure [3 top panel), and at the 
parallel and perpendicular temperatures (Figure [3 bot- 



tom panel, respectively in solid and dashed lines), we 
see that the anisotropy decreases straight away from the 
beginning, but not dramatically in value. The tempera- 
tures in Figure [7] are averaged over the whole active box. 
The decrease of the anisotropy is caused by the fact that 
the parallel temperature decreases faster than the per- 
pendicular one. At time TtUi ~ 40, however the parallel 
temperature reaches a plateau, causing an increase of the 
anisotropy, because the perpendicular temperature keeps 
decreasing. 

It interesting that the decrease of perpendicular tem- 
perature seems to be related to the geometry and the 
initial profile of magnetic field and density. Even if the 
initial mean velocity is zero, the magnetic field profile 
causes a narrowing of the pitch angle distribution for par- 
ticle moving outwards. Thus the initial conditions favor 
a distribution of particles that tend to align their velocity 
with the magnetic field, at the expense of the perpendic- 
ular velocity. From Figure [7] the parallel temperature 
will decrease if the firehose instability is triggered, but it 
might be that a more rapid decrease of the perpendicular 
temperature, when the linear stage of the instability has 
saturated, will result in a growth of anisotropy. Another 
interesting point is that the firehose instability, in this 
configuration, is not as effective in isotropize the par- 
ticles as is it would be for a non-cylindrical geometry, 
wi th constant magnetic field . Inde ed, it has been shown 
in \ Camporeale and Burgess] (|2008l ) , that in a Cartesian 
geometry, an anisotropic plasma is forced by the firehose 
instability to reduce its anisotropy to a state where the 
plasma remains close to marginal stability. 

4.3. Subsonic case: V m = 0.625 

In the subsonic case the interpretation of the results 
becomes less straightforward, because a vast propor- 
tion of electrons are now counter-propagating (i.e mov- 
ing towards the origin). This results in a non-locality of 
processes, where particles scattered at one location can 
rapidly influence the development of instabilities at other 
locations. Also, the convection of electromagnetic fluc- 
tuations towards outer regions is not as efficient as for 
the supersonic case. Figure |8] shows the development of 
5B/Bo, as for the previous cases. The whole dynamics is 
clearly slower, but the results are consistent with Figures 
[2] and [6j with again the formation of structures aligned 
with the background magnetic field. 

The top panel of Figure [9] shows the development in 
time of the temperature anisotropy (the format and leg- 
end are the same as in Figure These simulations are 
run again until Tlo.- l = 100, but the results should be 
compared with the results for the supersonic case, ob- 
tained only until Tuji ~ 25, since the flow speed is here 
four times slower. In the lower panel of Figure[9]we show 
with solid line the ratio of anisotropy for self-consistent 
run over test-particle run. The dashed lines show the 
same quantity, for the supersonic case (i.e. bottom panel 
of Figure [3]) , where now the time has been rescaled by a 
factor of 4. The two runs, for subsonic and supersonic 
flows, are qualitatively very similar, with the anisotropy 
reduced of about 5-15%. This is a sign that the growth 
rate of the firehose instability for an expanding plasma, 
must be dependent also on the flow speed of the plasma. 
Indeed, if this was not the case, the competition of the 
expansion and the instability would have led to quali- 
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tatively different results for the supersonic and subsonic 
cases. 

5. CONCLUSIONS AND DISCUSSION 

We have presented the results of PIC simulations of 
a plasma expanding radially on a disc, where the mag- 
netic field and density profile decrease linearly with the 
inverse of the distance. Although those simulations bear 
some simplification and assumptions with respect to a 
realistic stellar wind, they have the unique feature of 
treating self-consistently the effects due to the expan- 
sion and the electromagnetic fluctuations. In fact, we 
have used a fully-kinetic PIC code, with physical ion-to- 
electron mass ratio, and a computational box in a fixed 
frame of reference. Hence, we think that the results of 
this paper, summarised in this section, might be relevant 
for the understanding of a more realistic scenario. 

We have confirmed that the effect of electromagnetic 
fluctuations is to decrease the temperature anisotropy: 
while the expansion makes the parallel anisotropy grow, 
this increase is slowed down by the presence of EM fluctu- 
ations. It is interesting that the presence of fluctuations 
acts not only in decreasing the parallel temperature, as it 
is expected, but it also reduces the rate at which the per- 
pendicular temperature decreases during the expansion. 
The length and time scale constraints of our simulations 
limit our results to the evolution of the electron temper- 
ature anisotropy and its associated instabilities. 

The simulations have confirmed the presence of quasi- 
perpendicular waves, consistent with the development of 
firehose instability. However, the feedback effect played 
by fluctuations that counteracts the expansion, does not 
become more evident, when the plasma becomes linearly 
unstable, but is rather a continuous effect active since 
the beginning of the expansion. On the other hand, we 
have verified that if the expansion would not be present, 
and the plasma would be injected starting with a paral- 
lel anisotropy, this anisotropy would be reduced straight 
away. The fact that the results are not consistent with 
linear theory predictions should be thought as a conse- 
quence of both the geometry and the fact that the plasma 
is drifting. Both these effects are neglected in stan- 
dard linear theory, and although the approximation of a 



curved geometry with a planar one might be justified at 
large distances from the Sun, the drift should be proba- 
bly taken in account. A rough estimate of the importance 
of the drift, can be made by using the characteristic so - 
lar wind parameters listed 'm Wruno and Carbond {2005). 
The leading modes of the electron firehose instability 
have a wavevector that, depending on the anisotropy, 
ranges between kc/uji ~ 20 to kc/uji ~ 60. This corre- 
sponds, at 1 AU, to wavelengths of the order of about 10 
km. The growth rate is a function of the parallel beta, 
and the anisotropy, but if we take as a representative 
value of a fast growing mode a rate of 0.1S7 e (fi e is the 
electron cyclotron frequency) , then it would take approx- 
imately 5 ■ 10~ 2 seconds for the wave to grow by a factor 
of e. The bulk velocity of the solar wind varies between 
350 km/s (slow wind) and 600 km/s or more (fast wind), 
and the electron thermal velocity is between 2000 and 
3000 km/s. This means that the electrons can travel a 
distance comparable, if not greater, than the wavelength 
of the modes of interest, in a fraction of the growth time. 

Moreover, we have shown that the results for subsonic 
and supersonic flows are qualitatively similar, and the dy- 
namics of the subsonic flow is just slower. This suggests 
that the growth rate of the electron firehose instability 
must be a function of the drift speed of the plasma. This 
is because the expansion and the electromagnetic fluctu- 
ations play opposite role for the development of the tem- 
perature anisotropy. If the instability growth rate would 
not be a function of the drift speed, the temperature 
anisotropy would have been reduced more rapidly in the 
subsonic case, since here the increase of the anisotropy 
due to the expansion is slower. 

It has to be mentioned that another possible mecha- 
nism that is thought to isotropize the distribution func- 
tion is played by collisions. It has been reported that ex- 
ists an observational correlation between collisional age 
and electron temp erature anisotropy in the solar wind 
\Salem et all [2003) . Clearly, the role of collisions is not 
included in our simulations. 
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Fig. 2. — Supersonic case: 5B/Bo at times Tui = 18 (top), Tu>i = 60 (middle), Tu>i = 100 (bottom). The images are obtained by 
successive rotations of the box. 
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Fig. 3. — Supersonic case. Top panel: anisotropy T«/Tj_ at radial distances: ruii/c = 1.50 (red), ru)i/c = 1.86 (blue), ru^/c = 2.22 
(black). Solid lines are for the self-consistent run, and dashed lines are for the test-particle run. The dot-dashed line is the CGL prediction 
for ruii/c = 1.50. Bottom panel: ratio of anisotropy for self-consistent over test-particle runs. 




Fig. 4. — Supersonic case: Ty (top panel) and Tj_ (bottom panel). Solid lines are for the self-consistent run, and dashed lines are for the 
test-particle run, at radial distances: ru>i/c = 1.50 (r ed) , ru)j/c = 1.86 (blue), mij/c = 2.22 (black). The vertical dotted lines indicate the 
approximate time at which the condition of equation 1121 holds, that is the plasma is marginally stable. 



11 




Fig. 5. — Supersonic case. Electron distribution function at times Tui = 18 (top), To)j = 60 (middle), Tu)i = 100 (bottom), in (i)n,vj_) 
space. The distributions are averaged over the whole active box. 
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Fig. 6. — Static case: <5B/Bo at times Tu>i = 18 (top), Tu)i = 60 (middle), Tu>i = 100 (bottom). The images are obtained by successive 
rotations of the box. 
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Fig. 7. — Static case. Top panel: anisotropy Ty/Tx- Bottom panel: Ty (solid line), and T± (dashed line). 
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Fig. 8. — Subsonic case: 5B/Brj at times Tu)i = 18 (top), Tu)i = 60 (middle), Tu)i = 100 (bottom). The images are obtained by successive 
rotations of the box. 
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Fig. 9. — Subsonic case. Top panel: anisotropy T^/T± at radial distances: ru)i/c = 1.50 (red), ruii/c = 1.86 (blue), rui/c = 2.22 (black). 
Solid lines are for the self- consistent run, and dashed lines are for the test-particle run. Bottom panel: ratio of anisotropy for self-consistent 
over test-particle runs. Solid lines are for the subsonic case, and dot-dashed lines are for the supersonic case, with time rescaled by a factor 
of 4. 



